%% Main HF Organized

clearvars -GLOBAL
clear all
clc

rand('state',0);
randn('state',0);
% 1 - Tx.
% 2 - Channel Effects.
% 3 - Rx.
% ---------------------
% ---------------------
%% Initialization Parametrs
N_Pkts = 50;

Rate = 9600 ;                       % Supported rates are (3200-4800-6400-8000-9600) (Note: Rate 12800 is not supported yet)
Interleaver_L = 72;                % Interleaver sizes are (1-3-9-18-36-72)

ESNO_db = [20 24 28 33 35];

NOISE_FLAG = 1;                    % 0 for noise free and 1 for AWGN noise
Fractional_Flag = 1;                % zero for symbol space
SAMPLING_FLAG = 1;                 % set it to one for adding sampling offset effect
Phase_Offset = 1;
SAMPLING_Correction = 1;
N1_Synch = 25;
N2_Synch = 5;
Oversample = 4;

% Phase offset Upsampling ( used in solving the phase offset error).
PO_Upsampling = 16;             

% Carrier Frequency Offset, max frequency offset in the carrier is 0.018(Hz)
Freq_Offset = 0.0;

% The sampling offset shift between Tx and Rx (This will give 0.024 sampling error)
sampling_offset = -1/100000;
test_paramters = 0;
N1 = 45;                            % forward filter coefficients
N2 = 5;                              % backward filter coefficients
EsNo_test = (10.^(24/10));  % the for the SNR used in filter coefficients calculations (MMSE equation)


%% ************************************************************* %
% ******************* HF_INPUT_PARAMETERS *********************** %
% ************************************************************* %

Symbol_Rate = 2400;
Carrier_Frequency = 1800;
Sampling_Frequency = Symbol_Rate*Oversample;
SamplingFrequencykHz = Sampling_Frequency*1e-3;
% ************************************************************* %

% ************************** CHANNEL *************************** %
EsNo = 10.^(ESNO_db/10);

fadingMethod=1;
% fadingMethod:
% 0 = No fading
% 1 = Uniform Power
% 2 = exponentilal Power profile

pathDelaysMs = [0, 2];

% fadding (two segma) bandwidth
DopplerSpreadHz = 1.0;

% factor to compensate the equalizer transfer function at the high SNR
equalizerOffset = 0.0;


% for sampling error detection
Delay_Samples = 0;

% Number of symbols of the raised cosine filter.
Filter_Length = 12;

% Rolloff factor that satisfies the standard requirements on the filter (by trail)
Rolloff = 0.385;

% the number of the repeated Initial preamble (for AGC equirements)
% also we may not need to make AGC, but we have to detect
% the repeated sequence at the receive
N_Initial_Preamble = 1;

modulationIndex = find(Rate == [3200 4800 6400 8000 9600])+1;

% Constraint length of the encoder
K = 7;

% Code Rate of the encoder
codeRate = 3/4;
% ************************************************************* %

global LUT_PSK_Q_4 LUT_PSK_Q_8 LUT_QAM_16 LUT_QAM_32 LUT_QAM_64 Sync_Preamble Mini_Probe_Plus Mini_Probe_Minus LUT_S LUT_D LFSR_out Load_Location Mini_Probe_Plus_Modulated trellis_t LUT_Current
global interp_cons PFO_setimation;
PFO_setimation = [];
load Inputs/Preamble       % by (Preamble_LUTS.m), the output variables are (LUT_D,LUT_S,Mini_Probe_Minus,Mini_Probe_Plus,Sync_Preamble)
load Inputs/LUT_PSK_Q_4    % by (MODULATION_LUTS.m)
load Inputs/LUT_PSK_Q_8    % by (MODULATION_LUTS.m)
load Inputs/LUT_QAM_16     % by (MODULATION_LUTS.m)
load Inputs/LUT_QAM_32     % by (MODULATION_LUTS.m)
load Inputs/LUT_QAM_64     % by (MODULATION_LUTS.m)
load Inputs/alpha          % by (alpha_generation.m)
load Inputs/trellis_t      % by file (generate_trellis.m)
eval(['load Inputs/' 'LFSR_NSymbols_' num2str(modulationIndex)]);  % by (Gen_LFSRs.m), the variable output will be always (LFSR_out)
eval(['load Inputs/' 'Load_Location_Rate' num2str(Rate) '_Interleaver' num2str(Interleaver_L)]); % by (InterleaverLUTs_Generation.m), the variable output will be always (Load_Location)
LUT_Current = Get_Constellation(Rate);
load Inputs/Initial_Noise
load Inputs/init_noise_indecies
load Inputs/init_noise

disp('************************************************************')
disp('******************** INITIALIZATIONS ***********************')
disp('************************************************************')
disp(' ')
disp(['N_Pkts  =  ' num2str(N_Pkts)])
disp(['Rate  =  ' num2str(Rate)])
disp(['Interleaver_L  =  ' num2str(Interleaver_L)])
disp(['AWGN NOISE  =  ' num2str(NOISE_FLAG)])
disp(['ESNO db  =  ' num2str(ESNO_db)])
disp(['FADING  =  ' num2str(fadingMethod)])
disp(['DopplerSpreadHz  =  ' num2str(DopplerSpreadHz)])
disp(['Freq_Offset  =  ' num2str(Freq_Offset)])
disp(['sampling_offset  =  ' num2str(sampling_offset*SAMPLING_FLAG)])
disp('************************************************************')
disp(' ')

First_Sample_temp = [];
Paket_detec_temp =[];
%% Main Loop
for Pkt = 1:N_Pkts
    %% TX Transmitter Implementation
    % ########################################################### %
    % ########################################################### %
    % ################### Start Of TX ########################### %
    % ########################################################### %
    % ########################################################### %
    disp(['Packet Number = ', num2str(Pkt)]);
    disp('============================');
    %% INPUT DATA GENERATION
    % ----------------------
    IN_Frames = 72;
    Block_Size_In = HF_Get_Blocksize(Rate,Interleaver_L);
    Input_Bits = randi(1,[1 ceil(Block_Size_In*IN_Frames/Interleaver_L)]);
    % ----------------------------------------------------------- %
    % ----------------------------------------------------------- %
    %% BLOCKING
    % ---------
    [M_OUT_Blocks, N_Blocks] = HF_Blocking(Input_Bits,Rate,Interleaver_L);  % put the data into blocks
    Input_Bits_Actual=reshape(M_OUT_Blocks.',1,[]);
    % ----------------------------------------------------------- %
    % ----------------------------------------------------------- %
    %% BLOCK PROCESSING
    % -----------------
    Frame_Count = 1;
    for block_index = 1:N_Blocks
        OUT_Bits_Encoded(block_index,:) = HF_Convenc_Tailbiting_K7_Punct_R43(M_OUT_Blocks(block_index,:),trellis_t);%#ok
        OUT_Bits_Interleaved(Frame_Count:Frame_Count+Interleaver_L-1,:) = HF_Interleaving(OUT_Bits_Encoded(block_index,:),Rate,Interleaver_L);
        Frame_Count = Frame_Count+Interleaver_L;
    end
    N_Frames = Frame_Count-1;
    % ----------------------------------------------------------- %
    % ----------------------------------------------------------- %
    %% FRAME PROCESSING
    % -----------------
    for frame_index = 1:N_Frames
        OUT_Symbols(frame_index,:) = HF_Scrambling( OUT_Bits_Interleaved(frame_index,:),Rate);%#ok
        Modulated_Data(frame_index,:) = HF_Modulate(OUT_Symbols(frame_index,:),Rate);%#ok
    end
    % ----------------------------------------------------------- %
    % ----------------------------------------------------------- %
    %% FRAMMINGH
    % ----------
    TX_OUT = 1*HF_Framming(Modulated_Data,Rate,Interleaver_L,N_Frames,N_Initial_Preamble);
    % ----------------------------------------------------------- %
    % ----------------------------------------------------------- %
    %% Oversample
    % -----------
    TX_OUT_Oversampled = upsample(TX_OUT,Oversample);
    % ----------------------------------------------------------- %
    % ----------------------------------------------------------- %
    %% SRRC Filtering
    % --------------
    % Generate the Raised cosine filter
    HF_SRRC = rcosfir(Rolloff,Filter_Length/2,Oversample,1,'sqrt');
    TX_OUT_Filtered = conv(TX_OUT_Oversampled,HF_SRRC);
    Original_RX_IN = TX_OUT_Filtered;
    % ----------------------------------------------------------- %
    %%
    % ########################################################### %
    % ########################################################### %
    % ########################## End Of TX ###################### %
    % ########################################################### %
    % ########################################################### %

    %%
    % ########################################################### %
    % ########################################################### %
    % ################### Channel Effect ####################### %
    % ########################################################### %
    % ########################################################### %
    %% Fading Effect
    % ---------------
    [ RX_IN ] = WattersonFading( TX_OUT_Filtered, SamplingFrequencykHz, ...
        DopplerSpreadHz, pathDelaysMs, fadingMethod, Oversample, Fractional_Flag);
    % ----------------------------------------------------------- %
    % ----------------------------------------------------------- %
    %% Noise
    % ------
    noise = (randn(1,length(RX_IN))+1i*randn(1,length(RX_IN)))/sqrt(2);
    % ----------------------------------------------------------- %
    % ----------------------------------------------------------- %
    %%
    for SNR_Index =1:length(ESNO_db)
        tic
        disp(['SNR = ', num2str(ESNO_db(SNR_Index))]);
        disp('----------------------------------------------');
        
        %% Original Signal
        % -----------------
        Original_RX_IN = TX_OUT_Filtered;
        % ----------------------------------------------------------- %

        %% AWGN
        % ------
        Sigma = NOISE_FLAG / sqrt(EsNo(SNR_Index));
        RX_Faded_Noisy = RX_IN + Sigma/alpha(modulationIndex) * noise ;
        % ----------------------------------------------------------- %
        %% Sampling Frequency Offset
        % --------------------------
        if SAMPLING_FLAG == 1
            RX_Faded_Noisy = interp1(1:length(RX_Faded_Noisy), RX_Faded_Noisy,...
                1:(1+sampling_offset):length(RX_Faded_Noisy),'spline');
        end
        % ----------------------------------------------------------- %
        %% Phase offset effect
        % --------------------------
        if Phase_Offset == 1 
            Offset  = rand*sign(0.5-rand);
            RX_Faded_Noisy = interp1(1:length(RX_Faded_Noisy), RX_Faded_Noisy,...
                1-Offset:length(RX_Faded_Noisy),'spline');
            Offset_saved(Pkt, SNR_Index) = Offset;
        end
        % ----------------------------------------------------------- %
        %%
        % ########################################################### %
        % ########################################################### %
        % ################### Channel End ########################### %
        % ########################################################### %
        % ########################################################### %
        %% Rx
        % The Reciving Process
        % Time Synchronizatioin ,,
        
        % == >> Energy Detection ________________ Done
        % == >> Accurate Synchronization_________ Done
        % == >> SFO Correctioin & Estimation ____ Done with Good Results
        % == >> Sub-Carrir Mltiplication ________ Done
        % == >> Filtering With SRRC _____________ Done
        
        
        % Assuming No Phase Error Then Solving The SFO Error
        %% Energy Detection
        % ------------------
        [Paket_Start] = Energy_Detect(RX_Faded_Noisy);
        % ----------------------------------------------------------- %
        % ----------------------------------------------------------- %
        %% upsampling
        RX_Faded_Noisy = upsample(RX_Faded_Noisy, PO_Upsampling);
        Original_RX_IN = upsample(Original_RX_IN, PO_Upsampling);
        % ----------------------------------------------------------- %
        % ----------------------------------------------------------- %
        %% SRRC Filter
        % ------------
        OVS = PO_Upsampling * Oversample;
        HF_SRRC1 = rcosfir(Rolloff,Filter_Length/2,OVS,1,'sqrt');
%         HF_SRRC1 = HF_SRRC1/max(HF_SRRC1); % normlize filter
        HF_SRRC1 = HF_SRRC1*2^(log2(PO_Upsampling)/2);
        % Raised Cosine filter at the receiver
        RX_Faded_Noisy = conv(RX_Faded_Noisy,HF_SRRC1);
        
        % Original Tx Symbols passed through the Tx and Rx filters
        Original_RX_IN = conv(Original_RX_IN,HF_SRRC1);
        % ----------------------------------------------------------- %
        % ----------------------------------------------------------- %
        %% Time Synchronizatiton To Identify Start Symbol(course time synchronization)
        % ------------------------------------------------
        % Original_Signal = TX_OUT_Filtered;
        % The Original Signal At The Rx . End
        [First_Sample] = HF_Time_Synch(RX_Faded_Noisy,Original_RX_IN,...
            OVS,Paket_Start,HF_SRRC1);
        % Note :-
        %  We make  Rx*Carrier == >> Conv(Rx,SRRC)  == >> TO be Combatable
        %  With Previous System Only ,,, We Must Make Good Optimitation To
        %  this Function ,,, TO make It Lower at complixity
        % ----------------------------------------------------------- %
        First_Sample_Saved(Pkt, SNR_Index) = First_Sample;
                
        % ----------------------------------------------------------- %
        %% downsampling
        RX_Faded_Noisy = downsample(RX_Faded_Noisy(First_Sample+1:end), PO_Upsampling);
        Original_RX_IN = downsample(Original_RX_IN, PO_Upsampling);
        First_Sample = 0;
        % ----------------------------------------------------------- %
        %% SFO Correction And Estimation
        % -------------------------------
        % First Pkt We make Correction With Direct Estimation
        % After that all pkts will corrected thes update for next pkt the
        % estimation Value .....   
%         First_Sample = First_Sample + 24;
        Original_Start = 48;
        if (SAMPLING_Correction == 1)
            [RX_IN_Filtered] = SFO_Estimation(RX_Faded_Noisy, ...
                First_Sample,Oversample, N_Frames, Original_RX_IN,...
                Original_Start,Pkt);
        else
            RX_IN_Filtered = RX_Faded_Noisy ;
        end
        % ----------------------------------------------------------- %
        %% Information Extraction
        % ------------------------
        [Rate_Rx, Interleaver_Rx] = ...
            get_info(RX_IN_Filtered,Oversample, Original_RX_IN, EsNo_test, Original_Start,...
            N1_Synch, N2_Synch,First_Sample);
        % ----------------------------------------------------------- %
        % ----------------------------------------------------------- %
        %% Testing Variables for time_Synch & SFO_Corr &
        % Rate_Recived & Energy_detection & interleaver_length
        % ----------------------------------------------------------- %
        % ----------------------------------------------------------- %
        % ----------------------------------------------------------- %
        % ----------------------------------------------------------- %
        Paket_detec_temp (Pkt, SNR_Index) = Paket_Start;
        Rate_Saved(Pkt, SNR_Index) = Rate_Rx;
        Interleaver_Saved(Pkt, SNR_Index) = Interleaver_Rx;
        Actual_Start = init_noise(Pkt) + Original_Start;
        
        if(First_Sample == Actual_Start)
            FirstSample_Validate(Pkt, SNR_Index) = 1;
        else
            FirstSample_Validate(Pkt, SNR_Index) = 0;
        end
        
        if test_paramters == 1
            interp_cons
            continue ;
        end

        %% Equlization To Exrect Data
        % ----------------------------
        
        % First Iteration of Equlization
        % -------------------------------
            
        Start_Sample = First_Sample + (184 + 31 + 41) * Oversample;

        Hard = 1;                        % one to use hard decision, zero for all training
        [RX_Input_Frames, fadingSymbols, MSE_Real, MSE, OUT_Channel_OUT] = DFE_Equalizer(RX_IN_Filtered,...
            Start_Sample,Oversample,N_Frames,Rate_Rx,EsNo_test, Original_RX_IN,Hard,...
            N1,N2, Original_Start);
        
        [OUT_Bits_Decoded, DemodDescrambled_Data, N_Blocks] = HF_RECEIVER(RX_Input_Frames,N_Frames,...
            Rate_Rx,Interleaver_Rx,Sigma,fadingSymbols,alpha(modulationIndex),MSE_Real, MSE);
        
        Return_Tx_Data = reshape(OUT_Bits_Decoded,1,[]);
        
        [TX_OUT_Filtered,Input_Bits_Actual,OUT_Bits_Interleaved,N_Frames,TX_OUT] = HF_TRANSMITER(Return_Tx_Data,...
            Rate_Rx,Interleaver_Rx,N_Initial_Preamble,alpha,Oversample,HF_SRRC,modulationIndex);
        
        Original_RX_IN = conv(TX_OUT_Filtered,HF_SRRC);
            %% Second Iteration of Equlization
            % --------------------------------
        Hard = 0;            % one to use hard decision, zero for all training
        
        [RX_Input_Frames, fadingSymbols, MSE_Real, MSE, OUT_Channel_OUT] = DFE_Equalizer(RX_IN_Filtered, ...
            Start_Sample,Oversample,N_Frames,Rate_Rx,EsNo_test, Original_RX_IN, Hard, N1,N2, Original_Start);
        
        [OUT_Bits_Decoded, DemodDescrambled_Data, N_Blocks] = HF_RECEIVER(RX_Input_Frames,N_Frames,...
            Rate_Rx,Interleaver_Rx,Sigma,fadingSymbols,alpha(modulationIndex),MSE_Real, MSE);
        % ----------------------------------------------------------- %
        % ----------------------------------------------------------- %
        %% BER CALCULATION
        % -----------------
        % CODED BER
        % --------------
        RX_OUT_Bits = reshape(OUT_Bits_Decoded.',1,Block_Size_In*N_Blocks);
        biterrorvector = Input_Bits ~= RX_OUT_Bits(1:length(Input_Bits));
        BER(Pkt,SNR_Index) = biterr(Input_Bits,RX_OUT_Bits(1:length(Input_Bits)))/length(Input_Bits);
        
        AVG_BER(SNR_Index) = mean(BER(:,SNR_Index));
        
        BER
        AVG_BER
        
        toc
    end
    
end


save('Results\BER_Neq45_Ov16.mat', 'InterpConst_BER', 'AVG_BER');








